{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "298fcc05-0aff-4a5b-a1fb-1f1a04eda237",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "# 第一版"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2d5b4e45-b341-4ddd-a8a8-921639cc0cee",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ThermodynamicSystem:\n",
    "    def __init__(self, T=None, P=None, V=None, n=1, R=8.314):\n",
    "        \"\"\"\n",
    "        基础热力学系统类\n",
    "        T: 温度 (K)\n",
    "        P: 压强 (Pa)\n",
    "        V: 体积 (m³)\n",
    "        n: 物质的量 (mol)\n",
    "        R: 气体常数 (J/mol·K)\n",
    "        \"\"\"\n",
    "        self.T = T\n",
    "        self.P = P\n",
    "        self.V = V\n",
    "        self.n = n\n",
    "        self.R = R\n",
    "\n",
    "    def equation_of_state(self):\n",
    "        \"\"\"状态方程（需在子类中实现）\"\"\"\n",
    "        raise NotImplementedError(\"必须在子类中实现\")\n",
    "\n",
    "    def solve_variable(self, **kwargs):\n",
    "        \"\"\"求解缺失的热力学变量\"\"\"\n",
    "        for var, value in kwargs.items():\n",
    "            setattr(self, var, value)\n",
    "        return self.equation_of_state()\n",
    "        # 这个是从已知的值，然后返回一个方程？这样就用上sagemath解方程的本事了\n",
    "        # 后面看看子类里面怎么解的\n",
    "\n",
    "    def heat_capacity(self, process='constant_volume'):\n",
    "        \"\"\"热容（需在子类中实现）\"\"\"\n",
    "        raise NotImplementedError(\"必须在子类中实现\")\n",
    "\n",
    "    def entropy(self, T_initial, V_initial):\n",
    "        \"\"\"熵变计算（需在子类中实现）\"\"\"\n",
    "        raise NotImplementedError(\"必须在子类中实现\")\n",
    "\n",
    "    def enthalpy(self, T_initial, P_initial):\n",
    "        \"\"\"焓变计算（需在子类中实现）\"\"\"\n",
    "        raise NotImplementedError(\"必须在子类中实现\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3ff61e13-a12f-4d63-90d8-3de2d4e4fe7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "class IdealGas(ThermodynamicSystem): #这样写继承语法么\n",
    "    def equation_of_state(self):\n",
    "        \"\"\"理想气体状态方程 PV = nRT\"\"\"\n",
    "        vars_missing = sum(1 for var in [self.P, self.V, self.T] if var is None)\n",
    "        if vars_missing != 1:\n",
    "            raise ValueError(\"必须且只能缺失一个变量\")\n",
    "\n",
    "        if self.P is None:\n",
    "            return self.n*self.R*self.T / self.V\n",
    "        elif self.V is None:\n",
    "            return self.n*self.R*self.T / self.P\n",
    "        else:\n",
    "            return self.P*self.V / (self.n*self.R)\n",
    "\n",
    "    def heat_capacity(self, process='constant_volume', Cv=None):\n",
    "        \"\"\"\n",
    "        计算热容\n",
    "        Cv: 等容摩尔热容 (J/mol·K)\n",
    "        \"\"\"\n",
    "        if Cv is None:\n",
    "            # 默认单原子气体 Cv = (3/2)R\n",
    "            Cv = 3/2 * self.R\n",
    "        \n",
    "        if process == 'constant_volume':\n",
    "            return self.n * Cv\n",
    "        elif process == 'constant_pressure':\n",
    "            return self.n * (Cv + self.R)\n",
    "        else:\n",
    "            raise ValueError(\"无效的过程类型\") #我记得有个公式，把gamma作为变量的……这样的话绝热过程应该也能写\n",
    "\n",
    "    def entropy(self, T_initial, V_initial):\n",
    "        \"\"\"计算熵变 ΔS = nCv ln(T2/T1) + nR ln(V2/V1)\"\"\"\n",
    "        Cv = self.heat_capacity() / self.n  # 获取摩尔热容\n",
    "        delta_S = (self.n * Cv * ln(self.T/T_initial) + \n",
    "                   self.n * self.R * ln(self.V/V_initial))\n",
    "        return delta_S.n()\n",
    "\n",
    "    def enthalpy(self, T_initial, P_initial):\n",
    "        \"\"\"计算焓变 ΔH = nCp ΔT\"\"\"\n",
    "        Cp = self.heat_capacity(process='constant_pressure') / self.n\n",
    "        delta_H = self.n * Cp * (self.T - T_initial)\n",
    "        return delta_H.n()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "808d43bc-6f12-4347-9ab8-407b48b9d6ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "class RealGas(ThermodynamicSystem):\n",
    "    def __init__(self, a, b, **kwargs):\n",
    "        \"\"\"\n",
    "        范德瓦尔斯实际气体\n",
    "        a: 分子间吸引力参数 (Pa·m⁶/mol²)\n",
    "        b: 分子体积参数 (m³/mol)\n",
    "        \"\"\"\n",
    "        super().__init__(**kwargs)\n",
    "        self.a = a\n",
    "        self.b = b\n",
    "        #以sagemath的语法的话，应该会可以把a与b作为自变量来用的\n",
    "        # 而deepseek只是用python的语法了\n",
    "\n",
    "    def equation_of_state(self):\n",
    "        \"\"\"范德瓦尔斯方程 (P + a(n/V)^2)(V - nb) = nRT\"\"\"\n",
    "        vars_missing = sum(1 for var in [self.P, self.V, self.T] if var is None)\n",
    "        if vars_missing != 1:\n",
    "            raise ValueError(\"必须且只能缺失一个变量\")\n",
    "\n",
    "        if self.P is None:\n",
    "            return (self.n*self.R*self.T)/(self.V - self.n*self.b) - \\\n",
    "                   (self.a*(self.n**2))/(self.V**2)\n",
    "        elif self.V is None:\n",
    "            # 需要数值方法求解体积\n",
    "            raise NotImplementedError(\"体积求解需要数值方法\")\n",
    "        else:\n",
    "            return (self.P + (self.a*(self.n**2))/(self.V**2)) * \\\n",
    "                   (self.V - self.n*self.b) / (self.n*self.R)\n",
    "\n",
    "    def heat_capacity(self, process='constant_volume', Cv=None): # ERROR\n",
    "        \"\"\"实际气体热容（简化处理）\"\"\"\n",
    "        # 此处需要更复杂的模型，这里保持简单实现\n",
    "        if Cv is None:\n",
    "            Cv = 3/2 * self.R  # 近似值\n",
    "        return super().heat_capacity(process, Cv)\n",
    "\n",
    "    def entropy(self, T_initial, V_initial):\n",
    "        \"\"\"实际气体熵变计算（简化处理）\"\"\"\n",
    "        # 使用理想气体近似加上修正项\n",
    "        ideal_delta_S = IdealGas.entropy(self, T_initial, V_initial)\n",
    "        # 添加实际气体修正项（示例值）\n",
    "        correction = self.n * self.R * (self.b * self.n)/(self.V)\n",
    "        return (ideal_delta_S + correction).n()\n",
    "\n",
    "    def enthalpy(self, T_initial, P_initial):\n",
    "        \"\"\"实际气体焓变计算（简化处理）\"\"\"\n",
    "        # 使用理想气体近似加上修正项\n",
    "        ideal_delta_H = IdealGas.enthalpy(self, T_initial, P_initial)\n",
    "        correction = -2 * self.a * (self.n**2)/self.V\n",
    "        return (ideal_delta_H + correction).n()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4e2fa396-8e44-4813-b713-dc7c856624ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== 理想气体示例 ===\n",
      "压强: 199536.000000000 Pa\n",
      "等容热容: 24.9420000000000 J/K\n",
      "熵变: 5.43125117188022 J/K\n",
      "\n",
      "=== 实际气体示例 ===\n",
      "压强: 197885.558616704 Pa\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "ThermodynamicSystem.heat_capacity() takes from 1 to 2 positional arguments but 3 were given",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 14\u001b[0m\n\u001b[1;32m     12\u001b[0m real_gas\u001b[38;5;241m.\u001b[39mP \u001b[38;5;241m=\u001b[39m real_gas\u001b[38;5;241m.\u001b[39msolve_variable()\n\u001b[1;32m     13\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m压强: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mreal_gas\u001b[38;5;241m.\u001b[39mP\u001b[38;5;241m.\u001b[39mn()\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m Pa\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m---> 14\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m等容热容: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[43mreal_gas\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mheat_capacity\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mn()\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m J/K\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "Cell \u001b[0;32mIn[4], line 31\u001b[0m, in \u001b[0;36mRealGas.heat_capacity\u001b[0;34m(self, process, Cv)\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m Cv \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     30\u001b[0m     Cv \u001b[38;5;241m=\u001b[39m Integer(\u001b[38;5;241m3\u001b[39m)\u001b[38;5;241m/\u001b[39mInteger(\u001b[38;5;241m2\u001b[39m) \u001b[38;5;241m*\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mR  \u001b[38;5;66;03m# 近似值\u001b[39;00m\n\u001b[0;32m---> 31\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mheat_capacity\u001b[49m\u001b[43m(\u001b[49m\u001b[43mprocess\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mCv\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mTypeError\u001b[0m: ThermodynamicSystem.heat_capacity() takes from 1 to 2 positional arguments but 3 were given"
     ]
    }
   ],
   "source": [
    "# 示例用法\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=== 理想气体示例 ===\")\n",
    "    ideal = IdealGas(n=2, T=300, V=0.025)\n",
    "    ideal.P = ideal.solve_variable()\n",
    "    print(f\"压强: {ideal.P.n()} Pa\")\n",
    "    print(f\"等容热容: {ideal.heat_capacity().n()} J/K\")\n",
    "    print(f\"熵变: {ideal.entropy(280, 0.02).n()} J/K\")\n",
    "\n",
    "    print(\"\\n=== 实际气体示例 ===\")\n",
    "    real_gas = RealGas(a=0.365, b=4.28e-5, n=2, T=300, V=0.025)\n",
    "    real_gas.P = real_gas.solve_variable()\n",
    "    print(f\"压强: {real_gas.P.n()} Pa\")\n",
    "    print(f\"等容热容: {real_gas.heat_capacity().n()} J/K\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a92bbef2-a3c0-428c-8b08-35931854b8d8",
   "metadata": {},
   "source": [
    "# 第2版"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b751e66a-cd9e-4416-9392-22422be229c8",
   "metadata": {},
   "source": [
    "这版确实使用符号运算了\n",
    "不过还有点问题？（报的错我自己也不知道怎么解决）\n",
    "总之……试着拿来做题看看效果吧"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bb4b6781-ecad-4aa2-b1fa-4459655e2853",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sage.all import var, assume, solve, integrate, diff, ln, function\n",
    "\n",
    "class ThermodynamicSystemSymbolic:\n",
    "    def __init__(self, n=1, R=8.314):\n",
    "        \"\"\"\n",
    "        符号化热力学系统基类\n",
    "        n: 物质的量 (mol，符号或数值)\n",
    "        R: 气体常数 (J/mol·K)\n",
    "        \"\"\"\n",
    "        # 声明符号变量\n",
    "        self.T, self.P, self.V = var('T P V')\n",
    "        self.n = n\n",
    "        self.R = R\n",
    "        \n",
    "        # 假设所有物理量为正数\n",
    "        assume(T > 0); assume(P > 0); assume(V > 0)\n",
    "\n",
    "    def equation_of_state(self):\n",
    "        \"\"\"状态方程符号表达式（需子类实现）\"\"\"\n",
    "        raise NotImplementedError\n",
    "\n",
    "    def symbolic_relations(self):\n",
    "        \"\"\"显示基本热力学关系式\"\"\"\n",
    "        eq = self.equation_of_state()\n",
    "        Cv_expr = self.heat_capacity_symbolic('constant_volume')\n",
    "        Cp_expr = self.heat_capacity_symbolic('constant_pressure')\n",
    "        \n",
    "        print(\"状态方程:\")\n",
    "        show(eq)\n",
    "        print(\"\\n等容热容:\")\n",
    "        show(Cv_expr)\n",
    "        print(\"\\n等压热容:\")\n",
    "        show(Cp_expr)\n",
    "        return eq, Cv_expr, Cp_expr\n",
    "\n",
    "    def solve_symbolic(self, target_var):\n",
    "        \"\"\"符号求解指定变量\"\"\"\n",
    "        eq = self.equation_of_state()\n",
    "        return solve(eq, target_var, solution_dict=True)\n",
    "\n",
    "    def entropy_symbolic(self):\n",
    "        \"\"\"符号熵变计算 dS = (Cv/T)dT + (nR/V)dV\"\"\"\n",
    "        S = Function('S')(self.T, self.V)\n",
    "        dS_T = self.heat_capacity_symbolic('constant_volume') / self.T\n",
    "        dS_V = (self.n * self.R) / self.V\n",
    "        return dS_T*diff(self.T) + dS_V*diff(self.V)\n",
    "\n",
    "    def enthalpy_symbolic(self):\n",
    "        \"\"\"符号焓变计算 dH = Cp dT + (V - T(dV/dT)_P)dP\"\"\"\n",
    "        H = Function('H')(self.T, self.P)\n",
    "        dH_T = self.heat_capacity_symbolic('constant_pressure')\n",
    "        dVdT = diff(self.solve_symbolic(self.V)[0][self.V], self.T)\n",
    "        dH_P = self.V - self.T*dVdT\n",
    "        return dH_T*diff(self.T) + dH_P*diff(self.P)\n",
    "\n",
    "    def verify_maxwell_relation(self):\n",
    "        \"\"\"验证麦克斯韦关系 (dT/dV)_S = -(dP/dS)_V\"\"\"\n",
    "        S_expr = integrate(self.entropy_symbolic())\n",
    "        T_V = diff(S_expr, self.V).args[0].args[0]  # 提取温度项\n",
    "        P_S = diff(S_expr, self.T).args[0].args[0]  # 提取压强项\n",
    "        return bool(diff(T_V, self.V) == -diff(P_S, self.S))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b6acf4e6-2f06-4f78-ae87-a0cdf4912f9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "class IdealGasSymbolic(ThermodynamicSystemSymbolic):\n",
    "    def __init__(self, Cv_symbolic=False, **kwargs):\n",
    "        super().__init__(**kwargs)\n",
    "        self.Cv = var('Cv') if Cv_symbolic else (3/2)*self.R\n",
    "\n",
    "    def equation_of_state(self):\n",
    "        return self.P*self.V == self.n*self.R*self.T\n",
    "\n",
    "    def heat_capacity_symbolic(self, process='constant_volume'):\n",
    "        if process == 'constant_volume':\n",
    "            return self.n*self.Cv\n",
    "        elif process == 'constant_pressure':\n",
    "            return self.n*(self.Cv + self.R)\n",
    "        raise ValueError(\"无效过程类型\")\n",
    "\n",
    "    def entropy_symbolic(self, T_initial=None, V_initial=None):\n",
    "        \"\"\"ΔS = nCv ln(T2/T1) + nR ln(V2/V1)\"\"\"\n",
    "        if T_initial and V_initial:\n",
    "            return (self.n*self.Cv*ln(self.T/T_initial) + \n",
    "                    self.n*self.R*ln(self.V/V_initial))\n",
    "        return super().entropy_symbolic()\n",
    "\n",
    "    def isothermal_work(self):\n",
    "        \"\"\"等温过程功的符号表达式\"\"\"\n",
    "        return self.n*self.R*self.T*ln(self.V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "63f437d1-3f8d-4bab-8c52-257b1e3a7495",
   "metadata": {},
   "outputs": [],
   "source": [
    "class RealGasSymbolic(ThermodynamicSystemSymbolic):\n",
    "    def __init__(self, a, b, **kwargs):\n",
    "        super().__init__(**kwargs)\n",
    "        self.a, self.b = a, b  # 可以是符号或数值参数\n",
    "\n",
    "    def equation_of_state(self):\n",
    "        return (self.P + self.a*(self.n**2/self.V**2))*(self.V - self.n*self.b) == self.n*self.R*self.T\n",
    "\n",
    "    def virial_expansion(self, order=2):\n",
    "        \"\"\"维里展开式\"\"\"\n",
    "        #from sage.rings.puiseux_series_ring import PuiseuxSeriesRing\n",
    "        # 这是在sagemath里面，不用import也行吧\n",
    "        R = PuiseuxSeriesRing(QQ, 'V')\n",
    "        V = R.gen()\n",
    "        expr = (self.n*self.R*self.T)/(V - self.n*self.b) - self.a*(self.n**2)/V**2\n",
    "        return expr.truncate(order)\n",
    "\n",
    "    def heat_capacity_symbolic(self, process='constant_volume'):\n",
    "        \"\"\"使用状态方程计算热容差异\"\"\"\n",
    "        dPdT = diff(self.solve_symbolic(self.P)[0][self.P], self.T)\n",
    "        dPdV = diff(self.solve_symbolic(self.P)[0][self.P], self.V)\n",
    "        return self.n*self.R / (1 - (self.T*dPdT**2)/(self.P*dPdV))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e4b28708-74f0-450e-9f4a-46d1c09cd466",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== 符号运算演示 ===\n",
      "\n",
      "理想气体状态方程:\n",
      "状态方程:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html>\\(\\displaystyle P V = 8.31400000000000 \\, T n\\)</html>"
      ],
      "text/latex": [
       "$\\displaystyle P V = 8.31400000000000 \\, T n$"
      ],
      "text/plain": [
       "P*V == 8.31400000000000*T*n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "等容热容:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html>\\(\\displaystyle \\mathit{Cv} n\\)</html>"
      ],
      "text/latex": [
       "$\\displaystyle \\mathit{Cv} n$"
      ],
      "text/plain": [
       "Cv*n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "等压热容:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html>\\(\\displaystyle {\\left(\\mathit{Cv} + 8.31400000000000\\right)} n\\)</html>"
      ],
      "text/latex": [
       "$\\displaystyle {\\left(\\mathit{Cv} + 8.31400000000000\\right)} n$"
      ],
      "text/plain": [
       "(Cv + 8.31400000000000)*n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "熵变符号表达式:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html>\\(\\displaystyle \\mathit{Cv} n \\log\\left(\\frac{T}{T_{i}}\\right) + 8.31400000000000 \\, n \\log\\left(\\frac{V}{V_{i}}\\right)\\)</html>"
      ],
      "text/latex": [
       "$\\displaystyle \\mathit{Cv} n \\log\\left(\\frac{T}{T_{i}}\\right) + 8.31400000000000 \\, n \\log\\left(\\frac{V}{V_{i}}\\right)$"
      ],
      "text/plain": [
       "Cv*n*log(T/T_i) + 8.31400000000000*n*log(V/V_i)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== 实际气体分析 ===\n",
      "\n",
      "维里展开前两项:\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "unsupported operand parent(s) for -: 'Puiseux Series Ring in V over Rational Field' and 'Symbolic Ring'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[12], line 18\u001b[0m\n\u001b[1;32m     16\u001b[0m real \u001b[38;5;241m=\u001b[39m RealGasSymbolic(a\u001b[38;5;241m=\u001b[39ma, b\u001b[38;5;241m=\u001b[39mb, n\u001b[38;5;241m=\u001b[39mvar(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m'\u001b[39m))\n\u001b[1;32m     17\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m维里展开前两项:\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m---> 18\u001b[0m show(\u001b[43mreal\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvirial_expansion\u001b[49m\u001b[43m(\u001b[49m\u001b[43mInteger\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m     20\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m实际气体热容符号表达式:\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     21\u001b[0m show(real\u001b[38;5;241m.\u001b[39mheat_capacity_symbolic())\n",
      "Cell \u001b[0;32mIn[11], line 14\u001b[0m, in \u001b[0;36mRealGasSymbolic.virial_expansion\u001b[0;34m(self, order)\u001b[0m\n\u001b[1;32m     12\u001b[0m R \u001b[38;5;241m=\u001b[39m PuiseuxSeriesRing(QQ, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mV\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     13\u001b[0m V \u001b[38;5;241m=\u001b[39m R\u001b[38;5;241m.\u001b[39mgen()\n\u001b[0;32m---> 14\u001b[0m expr \u001b[38;5;241m=\u001b[39m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mR\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mT)\u001b[38;5;241m/\u001b[39m(\u001b[43mV\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mn\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mb\u001b[49m) \u001b[38;5;241m-\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39ma\u001b[38;5;241m*\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mInteger(\u001b[38;5;241m2\u001b[39m))\u001b[38;5;241m/\u001b[39mV\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mInteger(\u001b[38;5;241m2\u001b[39m)\n\u001b[1;32m     15\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m expr\u001b[38;5;241m.\u001b[39mtruncate(order)\n",
      "File \u001b[0;32m/usr/lib/python3.12/site-packages/sage/structure/element.pyx:1348\u001b[0m, in \u001b[0;36msage.structure.element.Element.__sub__ (build/cythonized/sage/structure/element.c:19623)\u001b[0;34m()\u001b[0m\n\u001b[1;32m   1346\u001b[0m     return (<Element>left)._sub_(right)\n\u001b[1;32m   1347\u001b[0m if BOTH_ARE_ELEMENT(cl):\n\u001b[0;32m-> 1348\u001b[0m     return coercion_model.bin_op(left, right, sub)\n\u001b[1;32m   1349\u001b[0m \n\u001b[1;32m   1350\u001b[0m try:\n",
      "File \u001b[0;32m/usr/lib/python3.12/site-packages/sage/structure/coerce.pyx:1282\u001b[0m, in \u001b[0;36msage.structure.coerce.CoercionModel.bin_op (build/cythonized/sage/structure/coerce.c:16906)\u001b[0;34m()\u001b[0m\n\u001b[1;32m   1280\u001b[0m     # We should really include the underlying error.\n\u001b[1;32m   1281\u001b[0m     # This causes so much headache.\n\u001b[0;32m-> 1282\u001b[0m     raise bin_op_exception(op, x, y)\n\u001b[1;32m   1283\u001b[0m \n\u001b[1;32m   1284\u001b[0m cpdef canonical_coercion(self, x, y):\n",
      "\u001b[0;31mTypeError\u001b[0m: unsupported operand parent(s) for -: 'Puiseux Series Ring in V over Rational Field' and 'Symbolic Ring'"
     ]
    }
   ],
   "source": [
    "# 示例使用\n",
    "if __name__ == \"__main__\":\n",
    "    print(\"=== 符号运算演示 ===\")\n",
    "    \n",
    "    # 创建理想气体符号系统\n",
    "    ideal = IdealGasSymbolic(n=var('n'), Cv_symbolic=True)\n",
    "    print(\"\\n理想气体状态方程:\")\n",
    "    ideal.symbolic_relations()\n",
    "    \n",
    "    print(\"\\n熵变符号表达式:\")\n",
    "    delta_S = ideal.entropy_symbolic(var('T_i'), var('V_i'))\n",
    "    show(delta_S)\n",
    "    \n",
    "    print(\"\\n=== 实际气体分析 ===\")\n",
    "    a, b = var('a b')\n",
    "    real = RealGasSymbolic(a=a, b=b, n=var('n'))\n",
    "    print(\"\\n维里展开前两项:\")\n",
    "    show(real.virial_expansion(2))\n",
    "    \n",
    "    print(\"\\n实际气体热容符号表达式:\")\n",
    "    show(real.heat_capacity_symbolic())\n",
    "    \n",
    "    print(\"\\n麦克斯韦关系验证:\", real.verify_maxwell_relation())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fea7fda7-bb3e-41e0-8cfe-5d3cc049769a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.5",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
